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Abstract. We show that the propagation of warps in 
gaseous disks can be strongly affected by compressional ef- 
fects, when the thickness of the disk is taken into account. 
The physical reason is that, in realistic self-gravitating 
disks, the sound time through the disk is comparable with 
the rotation time; thus the vertical hydrostatic equilibrium 
cannot be maintained adiabatically as the wave propa- 
gates (an implicit hypothesis in the thin-disk approxima- 
tion) and the disk cannot move up and down solidly to 
follow the warp perturbation. There results, together 
with the main vertical motion, a strong horizontal one 
which significantly modifies the dispersion relation. 

We then turn to the case of a disk composed of two 
fluids with different temperatures: this can correspond ei- 
ther to the combined motions of the gas and the stars 
in a galactic disk, or to their coupling with a flattened 
massive halo (assuming that, as for spiral waves, stars are 
conveniently represented by a fluid if one stays away from 
Lindblad resonances). We find, in addition to the usual 
warps, a short- wavelength wave which might explain the 
"corrugations" observed in many galactic disks. 

Finally, as a side result of this analysis, we discuss a 
possible weak amplification of m > 1 warps. 

Key words: galaxies: structure; spiral; kinematics and 
dynamics 



1. Introduction 

Warps are bending waves of disk galaxies, observed most 
often in the outer parts of edge-on spiral galaxies in the 
21 cm line of neutral hydrogen. Their physics is very sim- 
ilar to that of spiral density waves, from which they 
are distinguished by their opposite parity: they involve 
vertical rather than horizontal motions (though we shall 
see that horizontal motions in warps can be important as 
well) and antisymmetric (in the vertical coordinate) rather 
than symmetric perturbed potentials. Their propagation 
has been described in the limit of very thin disks, and 
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found to explain the main observational properties. Warps 
have also been considered as a possible source of angular 
momentum transfer in accretion disks. On the other hand, 
contrary to spirals, warps are not unstable, i.e. they can- 
not form spontaneously in an isolated disk. Since warps 
propagate radially in the disk, a continuous excitation 
mechanism must be found to explain their frequent oc- 
currence in galactic disks; various mechanisms have been 
considered (see e.g. Binney, 1992): 

— tidal excitation of warps (Hunter and Toomre 1969) by 
a companion galaxy; this cannot explain the observa- 
tion of isolated warped galaxies, and seems too weak to 
explain the observed amplitude of galactic warps. This 
has recently raised a strong interest in the context of 
accretion disks (Terquem 1993). 

— The bending of the trivial tilt mode by an elongated 
halo whose axis is misaligned with the axis of the 
galactic disk (Sparke 1984, see also Sparke and Caser- 
tano 1988). Although this hypothesis can describe very 
accurately the observed shape of prototype warps, it 
raises the problem of the origin and maintenance of 
such a misalignment. Recently Dubinski and Kuijken 
(1995) found that the axis of the halo and the disk 
should rapidly realign. 

This is the first in a series of papers where we will 
present a non-linear excitation mechanism for warps from 
spiral waves in disk galaxies. As a preliminary to this non- 
linear work we will concentrate here on the linear proper- 
ties of warps; we find important differences with the clas- 
sical theory, when we take into account the finite thickness 
of the disk. The reason is that the thin-disk approxima- 
tion used in the classical theory relies on an underlying 
hypothesis, that the whole gas column at a given location 
in the disk can move up or down solidly, without changing 
its vertical profile; this would require that the hydrostatic 
vertical equilibrium can evolve adiabatically, i.e. that the 
sound time through the disk, Ts = H/a where H is the 
disk scale height and a the sound speed, be shorter than 
the typical time scale of the warp, which is of the order 
of the rotation time tr = Jl"^. On the other hand the 
vertical hydrostatic equilibrium implies r^ ~ r^ if the ver- 
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tical gravity is dominated by a central object, or -a more 
appropriate condition for the outer galactic disk- if it is 
dominated by the self-gravity of the gas or the stars, pro- 
vided that the Toomre parameter Q is of the order of one. 

Finite thickness and compressibility effects have been 
intensively studied for spiral waves (Shu 1968, Vander- 
voort 1970, Romeo 1992, 1994). As an important result 
they found in particular that the critical Toomre parame- 
ter Qcrit for disk stability is lower than unity, as predicted 
by the thin disk analysis, at most by a factor about 2. 

Thus finite thickness and compressibility effects must 
also be taken into account for a realistic description of the 
warps. We find that they can substantially modify the dis- 
persion relation and that they introduce strong horizontal 
motions, comparable with the vertical one. 

Our method is similar to that of Shu (1968) for spiral 
waves, apart from modifications for taking into account 
the opposite parity of warps. 

On the other hand galactic disks are characterized by 
distinct populations, with different temperatures and scale 
heights, contributing separately to the vertical potential 
well: one has stars and gas, and also the halo whose contri- 
bution is believed to dominate the gravitational potential. 
In a first approach of the complexities introduced by these 
distinct populations, we present an analysis of warps in a 
disk composed of two fluids with distinct temperatures: as- 
suming that, as for spiral waves, stars can be conveniently 
represented as a fluid if one stays away from the Lindblad 
resonances, this allows us to discuss the case of warps in a 
disk of gas and stars, or in a disk embedded in a flattened 
halo. We find that, in addition to the usual warp wave, 
a new one exists which might correspond to the "corru- 
gations" observed in many galactic disks (Quiroga et al. 
1977, Florido et al. 1991). 

Finally we will also present an extension of this anal- 
ysis to a weak amplification mechanism for m > 1 warps. 
This is of very limited interest for galactic disks, where 
the amplification time would at best be a sizable fraction 
of the Hubble time, but might apply to accretion disks. 

2. Notations and formalism 

The dispersion relation of warps in an infinitely thin disk 
is (Hunter and Toomre 1969, Binney and Tremaine 1987) 



u^ 



fi^ + 27rGSg 



(1) 



where tu — a; — ?7i51(r) is the warp frequency in the frame of 
the matter, m is the azimuthal pattern number, /i is the 
vertical oscillation frequency of the disk, E is the mass per 
surface unit of the disk, and q the modulus of the warp's 
horizontal wavevector. 

We want to derive the dispersion relation for a moder- 
ately thick disk, i.e. a disk which is geometrically thin 
{H/r <C 1), where H is the disk thickness), but where the 
ratio ts/tt can be of the order of one. This will be done 
in a manner similar to the usual thin disk derivation, but 



we now have to consider the three projections (radial, az- 
imuthal and vertical) of the Euler equation, and the solu- 
tion of the Poisson equation becomes more complex. 

2.1. Notations 

Let us first introduce some notations. We use the shearing 
sheet model to describe differential rotation (Goldreich 
and Lynden-Bell 1965). In this model one considers an 
annulus around corotation as a cartesian slab with x = r — 
ro (where vq is the corotation radius where a; = 0) and y — 
r-ff. The radial variations of all equilibrium quantities are 
neglected, except the rotation speed Vo{x)ey which varies 
linearly, Vq = ro^{ro) + 2Ax where A = l/2rdft/dr is 
Oort's first constant. We also restrict ourselves to the case 
of isothermal equilibrium and perturbations, with sound 
speed a. 

In what follows we denote by U, V and W the com- 
ponents of the perturbed velocity, p the perturbed den- 
sity and (j) the perturbed potential. We call kx and ky 
the projections of the wavevector of the warp we study, 
Q = i^x + ^yY^'^ its modulus, Co its frequency. The con- 
nection between the cylindrical geometry and the shearing 
sheet is given by ky = m/rQ. Unperturbed quantities are 
denoted by the same symbols with a subscript 0. 

2.2. The formalism 

The dispersion relation is derived in a WKB approxima- 
tion, whose validity condition is the "tightly wound" ap- 
proximation, kx ^ ky. In section 5 we will give a numerical 
solution, independent of this approximation. 

We limit ourselves to a linear analysis of the propaga- 
tion of warp waves, i.e. perturbed quantities are infinites- 
imal. With these restrictions, it is possible to characterize 
a warp wave by parity considerations, since the parity of 
the equilibrium and of the perturbation equations allows 
one to separate symmetric and antisymmetric solutions. 
In the former, corresponding to spiral density waves, all 
perturbed quantities are even in z, except W which is 
odd. The opposite holds for antisymmetric solutions, cor- 
responding to warps. 

2.3. The basic equations 

We start from the continuity equation and the horizontal 
projections of the Euler equation: 



d_ 



iLop + po{ikxU + ikyV) + —^{pnW)^Q (2) 



iu)U — 2VtV ~ —ikx \ (f) + a 
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where k = (417^ + 4riA)^/^ is the epicychc frequency. 
Using the WKB assumption (kx ^ ky) we get: 



9 UJ 

kxU + kyV = q-'^TT, 



Po 



Thus the continuity equation (g) can be written as: 



da 
dz 



where 

a = ipo{z)W 

P 

s = — 

Po 



1 ~ 
q^'uj 



' + 0^8) 



Poiz) 



(5) 



Our basic set of equations also includes the vertical 
projection of the Euler equation: 



ds _ 1 
dz a^ 



bja d<p 
Pa dz 



(6) 



and the Poisson equation: 



dz 



-^-^ = 47rG'/9os + q 



In equation (^ the second term in the bracket rep- 
resents the horizontal part of the divergence of the ve- 
locity field, i.e. the compressional contribution which is 
the main new physical effect introduced in this paper. For 
J) o^ 51 ^ Kj and unless (D is very close to k {i.e. at the 
Lindblad resonances), this term is of order oiq^a^ / k^\ this 
is of the order of q^H^ compared to the first term, if the 
disk is at hydrostatic equilibrium, with a gravity either 
dominated by a central object or (if Toomre's parameter 
Q — na/irGYi ^ 1) by the local gravity of the gas, so that 
H ^ a/ K. We will return to this point in the following sub- 
section. Restricting ourselves to wavelengths larger than 
the disk thickness, we will use qH as the expansion pa- 
rameter in the following analysis. 

On the other hand, we will ignore here the case, which 
might apply to galaxies, where gravity is dominated by a 
massive halo, so that H ^ a/ k. If the halo is passive, i.e. 
does not participate in the motion, the disk is indeed thin 
in the sense that ts ^ tr. On the other hand, if the halo 
is flattened it presumably participates in the differential 
rotation and can also be involved dynamically in the warp. 
Its effect will be discussed below, in the section devoted 
to warps in a two-component disk. Even in the case of 
a passive halo, we will find new effects associated with 
corrugations of the disk. 



2.4. Consistent model oj geometrically thin disk 

In the following we will make extensive use of consistent 
models of disk vertical density and potential profiles. We 
build them as coupled solutions of the hydrostatic equi- 
librium and the Poisson equations, in the hypothesis of a 
geometrically thin disk [H/r <C 1). 

More precisely, let us construct ah nihilo such a con- 
sistent disk. In a first guess, we choose the density profile 
pt){z). The vertical hydrostatic equilibrium gives the grav- 
itational potential by: 



Po- 



dz 



2 9 , . 
= a -poiz) 



On the other hand the Poisson equation gives: 



1 d 



d 



dr V dr 



h = 4:TtGpo - 



d^o 
dz"^ 



where we have written what we know [i.e. what is 
imposed by our choice of po) in the R.H.S. Thus the Pois- 
son equation gives the behavior of the L.H.S., which is the 
radial part of the laplacian, hereafter denoted by A^: 



1 d 



A,. = -^ r— 



d 



r dr \ dr 



It is an easy matter to check that this radial laplacian 
is equal to: 



(7) A^(/.o = K^ -2V? = - 



P 



(8) 



Note that /i^, that we have called the vertical oscil- 
lations frequency, does not involve any term in '^■KGpm. 
Indeed AnOpm appears when considering the vertical os- 
cillations of a test particle in the rest potential of the 
disk, whereas we are concerned here with global oscilla- 
tions which involve vertical motion of the potential well 
itself. (see Hunter and Toomre, 1969, for a derivation 
of equation S). Since the density profile is arbitrary, so 
is p?' . But in a geometrically thin disk Vl and k must be 
independent of z. Thus A^^o must also be constant with 
respect to z. We can thus summarize the definition of a 
consistent disk equilibrium as a disk whose density profile 
obeys the equation 



,d_ 
dz 



1 dpf) 
Po dz 



AnGpo 



(9) 



Where both a and Aj.0o are independent of z. 

By choosing them we can construct a continuous 
palette of vertical equilibria, from a keplerian disk (where 
the vertical gravity is dominated by the mass of a cen- 
tral object, i.e. where the first term is dominant in the 
R.H.S. of equation M) to self-gravitating ones, where 
gravity is dominated by the local distribution of mat- 
ter {i.e. where the second term is dominant in equation 
yi). One easily obtains in the first case a gaussian profile 
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Poiz) = pmexp{—z'^/H'^) where jjl^ = —Ar4>o = 2a^/_ff^ 
(and in that case, /i and k are equal to fi). The second 
case straightforwardly leads to: poiz) — pm. / cosh^{z/ H) 
where H — a/y/2nGp^. Between these extremes, one can 
consider any mixture of local and global gravity. 

Since the construction of a consistent disk implies 
the choice of two independent parameters (the sound 
speed a and p,"^), the set of consistent disk equihbria is 
a two-dimensional manifold where any independent cou- 
ple of parameters can be chosen to "label" a particular 
choice; We will use below p/n and the Toomre parameter 
Q — aK/nGT,. 

Then obviously the notion of consistent disk equilib- 
rium does not imply any constraint on the Toomre param- 
eter Q, and on the ratio ts/tr discussed in the introduc- 
tion. 

However, in order to build realistic disk equilibria, one 
cannot choose p/ k and Q arbitrarily: indeed if p/ k is close 
to one, meaning that the disk is nearly keplerian and dom- 
inated by the gravity of the central regions, the Q param- 
eter must be large enough to ensure that the local gravity 
is low enough. On the other hand, if p/ k is not close to 
one, Q must not be too large, so as to enable the local 
gravity to play a role in the equilibrium of the disk. Let 
us also recall, as already mentioned, that we will not dis- 
cuss here the contribution of a halo to the vertical gravity 
and thus to the frequency p. This discussion is deferred to 
the consideration of two-component disks in section 4. 

In the following we study the dispersion relation of 
warps for various type of disks, from self-gravitating to 
keplerian. We label them by the quantity p/ k. Of course 
we have checked that the Toomre parameter of these disks 
is realistic in the sense that it is low for self-gravitating 
disks and very large in the case of the keplerian disk. 



A fifth condition corresponds to the amplitude of the so- 
lution (arbitrary since we consider linearized equations). 
Since we solve a fourth-order system, these conditions can 
be fulfilled only if the parameters q and uj obey a con- 
dition, which is the dispersion relation. This dispersion 
relation can also be seen as a compatibility condition of an 
overdetermined differential system which has both Dirich- 
let's and Neumann's boundary conditions, as explained in 
Bertin and Casertano (1982), who investigated the dis- 
persion relation of bending waves in incompressible thick 
disks with constant thickness. 

3.1. Numerical determination of the dispersion relation 

We determine numerically the dispersion relation by solv- 
ing the system (5-7) and looking for the values of q and 
Cj which allow the solution to obey the boundary condi- 
tions. Our procedure is the following: for a given value of q 
and choosing an approximate value of d) we integrate the 
system (5-7) from large z to zero. Obeying the bound- 
ary conditions at infinity means that we start with a = 
and d(j)/dz = —90 (note that an error in the latter con- 
dition will give a projection on the solution that diverges 
at large z, so that as we integrate toward decreasing z its 
effect will be exponentially small); we choose s — 1 and we 
still have one free parameter, the value of 0. We do this 
twice, starting with different values of 0, and resulting in 
two different sets of values of s and (f) aX z = 0. We can 
then find a linear combination of these two solutions to 
form a third one, which still obeys the conditions at in- 
finity and now also has s(0) = 0. Then a Newton method 
allows us to vary the value of w, for a given g, until the 
last condition, 0(0) = is also satisfied. This results in 
Co{q), the dispersion relation. 



3. The dispersion relation of v^rarps 

We write the dispersion relation by solving the linear sys- 
tem (5-7) for a, s and 0, subject to the appropriate bound- 
ary conditions. We have four such conditions, two at the 
mid plane ensuring the parity of the solution: 

s = 

= 

and two at infinite z: 

a = 

oz 

which ensure that the mathematical solution which 
satisfies the warp parity is also a physical one, with an ex- 
ponential decrease of the perturbed potential, and a van- 
ishing momentum fiux due to the rarefaction of matter. 



3.2. Analytical derivation of the dispersion relation 

We remind that equation (^ has been obtained for an in- 
finitely thin sheet of matter involving only vertical motion. 
Thus this dispersion relation does not take into account: 

— compressional effects due to the finite sound speed 
(the finite ratio ts/tb). We will emphasize this point 
further below, when we derive the analytical dispersion 
relation for a moderately thick disk; 

— horizontal motions; 

— geometric effects, often approximated by "softened 
gravity" models for spiral density waves. Analytical 
calculations for an incompressible disk in which hor- 
izontal motion is suppressed show - for any density 
profile - that Cj tends to an asymptotic value Coao when 
qH becomes large which is obviously not the case with 
equation ([^). 

We present a method which allows us to derive ana- 
lytically an approximate dispersion relation, taking into 
account the compressional effects and the finite thickness 



F. Masset & M. Tagger: Propagation of warps in moderately thick disks 



for small but finite values of qH. It must be emphasized 2 
that the thin-disk dispersion relation gives . 

for (D ~ ri, and a disk which has either a low self- 
gravity, or a strong one with the Toomre parameter Q = 
Ka/nGT, ^ 1 (the latter case being more relevant for the 
outer regions of disk galaxies). Thus qH ^ QH/a — Ts/tu, 
the ratio of the sound time through the disk to the dynam- 
ical time: the thin disk approximation (H —f 0) involves 
an assumption that the sound time through the disk is 
small, so that the vertical hydrostatic equilibrium can be 
maintained throughout the evolution of the waves of in- 
terest and compressional effects can be neglected. On the 
contrary the new effects we find here, which are associated 
with the finite value of qH, result from the fact that this 
equilibrium cannot be maintained perfectly. 

Our method, which aims at eliminating the z depen- 
dence of the variables, consists in making an integral over 
z which represents the amplitude of the warp (the mean 
displacement from the galactic plane) , and transforming it 
by using equations (0-0) • The main difficulty comes from 
the potential. Here we express it by using Green's func- 
tions, which allows us to separate the vertical structure of 
the solution without having to average over z, as is clas- 
sically done. We show that it is possible to modify the 
dispersion relation by a second order term in qH. This 
term, similar to the pressure {a?q^) term of the disper- 
sion relation of spiral modes, contains the contribution of 
the compressional effects, but does not contain the mod- 
ification of the gravity force, which is more difficult to 
evaluate, and would require the analytical knowledge of 
the eigenvectors. 

The derivation can be found in appendix, and it leads 
to: 




h ig. 1. In this figure we have assumed a purely solid vertical motion. 
Hence the density is conserved and advected by the vertical motion, and, 
since pressure is linked to density by: P — a p, horizontal gradients of 
density and pressure arise in the disk, due to the vertical stratification 
of the equilibrium configuration. These gradients are proportional to q, 
the horizontal wavenumbcr, and to H~ , the vertical gradient. 



u' 



= -Ar 



27rGS9 - 



2 2 
ja q 



Ograv[{qHf 



are usually very similar, except in the vicinity of the Lind- 
blad resonances where a major difference occurs between 
the two cases: the stars resonate with the wave and can ex- 
change energy with it through Landau damping, whereas 
the gas does not. Here the divergence of this term does not 
mean a resonance, but rather the fact that our expansion 
fails in this region. Still, we note that this allows us to 
get close enough to the Lindblad resonance that our term 
starts playing an important role, for stars as well as for 
the gas. 

We also wish to note that a similar contribution to 
the dispersion relation had been found by Nelson (1976 
and 1981) in the limit of weak shear, and more recently 
by Papaloizou and Lin (1995) in a more general study in 
cylindrical geometry. 



where Ograv[{qHY] represents the (unknown) contri- 
bution of self-gravity at orders 2 and higher in qH . 

Using equation (|8|), and neglecting the second-order 
gravitational contribution, we finally obtain: 



~2 2 

U = n 



2TiGY.q 



-,a\\ 



(10) 



This dispersion relation is identical to the thin-disk 
one (^, except for the second-order compressional term. 
This term comes from the horizontal motions driven by 
the horizontal gradients of the perturbed pressure and 
potential, and diverges at the Lindblad resonance where 
epicyclic motion is resonantly excited by the pressure 
force. However this divergence must be taken with cau- 
tion: we have thus far discussed only the case of a gaseous 
disk, through the use of hydrodynamical equations. It is 
well known that the results for a coUisionless stellar disk 



3.3. Physical interpretation of the additional term 

The additional term in the dispersion relation is linked 
to the divergence of the horizontal perturbed velocity (as 
can be seen in the derivation given in appendix, from equa- 
tion [^). Hence this additional term is due to the pres- 
ence of horizontal, compressional motions associated with 
the warp. In order to understand how these horizontal 
motions arise, let us first consider a warp in which motion 
is purely vertical and solid. This situation is presented in 
figure H 

It shows how horizontal pressure gradients appear in 
the system, which in turn tend to move matter horizon- 
tally. The amplitude of horizontal motions is greater and 
greater as the warp frequency in the frame of matter (w) 
approaches the frequency at which matter spontaneously 
moves horizontally, i.e. the epicyclic frequency. This can 
be the beginning of a justification of the resonant denom- 
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r" ig. 2. Wg show (left) the velocity field (U^W) of a warp far from 
the Lindblad resonance, and (right) the same field close to the Lindblad 
resonance. If we define R as the ratio of the 27rGS(y term in the dispersion 
relation over the additional a q uj /{uj — n ) term, we have ii — 41 in 
the first case, and R — 0.53 in the second case. These fields have been 
obtained numerically by the method described in the text. The disk scale 
height H is about 2 in these plots. 



inator in the additional term to the dispersion relation. 
In order to illustrate the behavior of matter far and close 
to the Lindblad resonances, we have plotted in figure H 
the velocity fields of matter in these two cases. We see 
that far from the Lindblad resonance, the matter moves 
essentially vertically with the motion described by the dis- 
persion relation without our additional term; close to the 
Lindblad resonance, on the other hand, the motion tends 
to be practically horizontal as soon as we leave the me- 
dian plane. In this case, the motion might be described by 
two slices moving in opposite directions, giving a strong 
vertical shear. 

In a forthcoming paper discussing the possible excita- 
tion of warps by spiral waves, we will fully analyze the 
energetics of warps, and in particular we will discuss the 
fraction of the total energy of the warp that can be stored 
in horizontal motions. 

It is noteworthy that our additional term is linked to 
the finite thickness of the disk, since a fundamental hy- 
pothesis is that the disk equilibrium is consistent, and 
since the sound speed which appears in that additional 
term is proportional to the thickness of the disk. 



3.4- Comparison between analytical derivation and 
numerical results 

We present the comparison between the numerical disper- 
sion relation and the dispersion relation given by equa- 



tion (Jl^) in figure ^. We obtain a very good agreement 
between the numerical and analytical solutions, which ap- 
pears to be far better than the infinitely thin disk disper- 
sion relation (m, except in the top- left diagram, where the 
disk is strongly self-gravitating, so that the gravitational 
part of the second-order term, which we did not evaluate, 
cannot be neglected. 

We have also emphasized in this figure the case of a 
massless keplerian disk which is for us, according to our 
definition of consistent disk, a massless disk (S = 0), for 
which il ex r~^'^. For such a disk, the general dispersion 
relation becomes: 



to' 



K ± aujq 



This dispersion relation is exact , since the unknown 
correction term in our derivation was due to self-gravity, 
which we do not need to take into account for a keplerian 
disk. Thus we can see that the warp of a keplerian disk 
precesses, a result which was not given by the infinitely 
thin disk relation which in that case becomes: 



c:;2 






The precession frequency, expanded to lowest order in 
qH, is thus of the order of 



riK 



H 

'r 



(for a wavelength of the warp of the order of the size of 
the disk). This can have important consequences in the 
context of protoplanetary disks. The role of warps has 
recently raised a strong interest, based in part on the fact 
that in keplerian disks the warp is a neutral tilt mode 
which has a vanishing energy and can thus very easily be 
excited. We expect that our effect should involve a finite 
warp energy, and thus make excitation mechanisms less 
efficient. 

A similar resolution for the spiral waves dispersion re- 
lation (where we omit again the second and higher orders 
terms in qH for the self-gravity) leads to the infinitely thin 
disk dispersion relation; 



= K^ - 27rGSq 



2 2 

a q 



(11) 



From this result we can say that equation (hlf) and 
equation (|l^ describe respectively the propagation of spi- 
ral waves and of bending waves with the same accuracy, 
i.e. with the same underlying hypothesis. 

4. Tvifo-fluid dispersion relation 

In this section we derive the dispersion relation of the 
bending waves in a two-fluid system. Both fluids have dif- 
ferent sound speeds and hence scale heights, and they have 
independent surface densities. Our primary goal is to dis- 
cuss warps in a disk composed of stars and gas (assuming 
that the stars are correctly represented as a warm fluid). 



F. Masset & M. Tagger: Propagation of warps in moderately thick disks 



where we have defined: 




0.01 0.02 0.03 0.04 0.05 
qH 



qH 



r" ig. 3. Comparison between the analytical and numerical dispersion 
relations. The quantity v represents Cj j ti. The triangles represent the 
results of numerical integration, the dashed line the infinitely thin disk 
dispersion relation, and the solid line the thick disk dispersion relation. 
For each curve we have chosen a range of qH such that the diverging 
trend between the numerical dispersion relation and the infinitely thin 
disk dispersion relation appears at about the first quarter of the qH 
range. The ordinate in the bottom right diagram is v rather than v . 
See text for explanations. 



or to a disk embedded in a flattened halo which can par- 
ticipate dynamically in the warp. In the following, for sim- 
plicity, we refer to the two fluids as gas and stars. 

We will find in particular that, in a star-gas disk, 
we find solutions closely resembling the corrugations, i.e. 
short wavelength bending waves of the gas layer observed 
in certain edge on spiral galaxies and in the Galaxy. 

We keep the same notations as for the mono-fluid case, 
with an index * or g applying respectively to the warmer 
(stars) and the cooler (gas) species. 

We still need the assumption of consistent vertical 
equilibrium already used in our mono- fluid derivation, i.e. 
that both fluids obey the Poisson equation and the hy- 
drostatic equilibrium, while the gravitational potential is 
such that Ar^o is constant through the disk thickness. 

With a method strictly similar to that exposed in ap- 
pendix for the monofluid case, we can write down equa- 
tions for Nf and Ng, which constitute a linear homoge- 
neous system. This system will have a non-trivial solution 
only if it has a vanishing determinant. 

Hence the dispersion relation of warps in a two-fluid 
system is: 



VgV, - . 



/27rGE„g 



/27rGS,jg 



= 



(12) 



1 m' 
Si 



27rGSig 



where i — * or g 



and: 



= AttG 



r + oc 
Jo PjO'^ 

N,, 



where j y^ i 



Physically, Vi = would be the dispersion relation of 
bending waves if there was only the species i. The role of 
I'g becomes clear if one goes to the limit where the gas 
disk is much thinner than the stellar one: in that case one 
has Vg = AirGprn*, where pm* is the stellar density at the 
disk mid-plane; Vg then appears as the vertical oscillation 
frequency of the gas in the stellar potential. 

It must be noted that realistic stellar disks can have a 
vertical velocity dispersion very different from the radial 
one. This might affect the present results by numerical 
factors which could be important in detailed comparisons 
with observations. This will be considered in future work. 

For a given (D, equation (^2|) is of order 4 in qH (since 
the 1/Si hidden in the Di are of order 2 in qH). In general 
its roots cannot be easily expressed, but we find it inter- 
esting to first analyze them in the trivial case ag — a,, 
before turning to numerical solution in the general case. 

4.1. Roots identification of the two-fluids dispersion 
relation 



4.1.1 



Peculiar case ag — a^. 



In this case we have 6*3 = 5* = S, and we are dealing 
with two physically indistinguishable fluids. Now we can 
factorize the dispersion relation (|lj) . We obtain after some 
straightforward transformations: 



^2 + 27rGStg 



where Et = S.^ + Eg. The eigenvector associated with 
the second factor can be easily found by summing equa- 
tions relative to iV, and Ng, giving: 



^2 + 27rGEtg 



{Ng + N,)=0 



Thus when the second factor vanishes this gives: Ng + 
N^, = 0. This condition consists in splitting the gas in two 
parts that move in opposite directions, so that the "total" 
warp (the average displacement) vanishes. In fact we note 
that the second factor is of the form aq'^ + b, so that it gives 
only one positive root for q. We will call this the "hidden" 
mode, even when we relax the condition ag ~ a^. 

In the same manner, it is an easy matter to check that 
the eigenvector of the flrst factor is E^A'g — E^A^, — 0. 
This means that both fluids have the same behavior (they 
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r ig. 4. The diniensionless wavenuinber of the hidden mode for values 
of parameters typical of our Galaxy's, and the ratio Ng/N^, which is 
ten times (S*/Sg) the ratio of deviations of the stars and gas. The ratio 
N^/Ng vanishes at ag/a» — 0.6524 (i.e. at qH^ — 2), which corresponds 
to a regime where the stars are strictly motionless. 



have in particular the same elongation Z) , and we recover 
the one-fluid dispersion relation. 

Finally we have three modes: the two "classical" modes 
of the one fluid disk, given by the second order dispersion 
relation of a one-fluid warp, and the "hidden mode" asso- 
ciated to the second factor. 

4.1.2. General case 

When Qg and a* are different, the simple factorization 
seen above is no more possible and we turn to numerical 
solution. 

Wc have adopted the following values: 

fi^/n^ = 0.1 

(corresponding to a nearly flat rotation curve) 

1/2/^2 ^ 10 

Cj/k = 0.75 

(so that we are between the warp's Lindblad resonance 
(inner or outer) and the forbidden band around corotation 
where the warps do not propagate) 

Eg/E* =0.1 

We also need to choose the value of i/*. It is an easy 
matter to check that, since we limit ourselves to the (re- 
alistic) case Og ^ a*: 



r ig. 5. Plot of the relative error between the exact numerical solution 
and the approximate dispersion relation. We see that the error remains 
weak, and vanishes when agja^, — 0.6524, which is logical since at this 
point the stars are strictly motionless. 



Vg ~ 47rGE*/iJ, and i/* ~ AnCEg/H^ hence z^» ~ 

We will not discuss the standard warp mode, which 
behaves in a manner very similar to the one-fluid case. The 
stars and gas motions are nearly identical, though since 
flg j^ a* small differences occur. We will rather discuss the 
behavior of the "hidden" mode, and tentatively identify it 
with the corrugations observed in many galaxies. 

The results are plotted on figure ^. The top plot shows 
qH^,, while the bottom one shows the ratio N^,/Ng, as a 
function of Ug/a^, for the hidden mode. When the sound 
speeds are identical we check the result of the previous 
subsection, that the two fluids move in opposite direc- 
tions so as to achieve a null global vertical displacement. 
When the gas sound speed is decreased, the ratio |A''*/A'^g| 
decreases, i.e. the stars are more and more motionless. 
Thus the mode under study is essentially gaseous. The 
stars become passive, though they still act through their 
unperturbed potential to change the frequency of verti- 
cal oscillations of the gas disk. One may note that for 
ttg/a^, < .65 the ratio N^/Ng changes sign, i.e. the stars 
now move in the same direction as the gas. 

We have monitored the relative error obtained with the 
computation of q by the approximate dispersion relation: 



Cj^ 



iiiGYsgq 



2 2 



i.e. the one-fluid dispersion relation where we have 
added the vertical frequency Vg due to the rest potential 
of the stars. 



F. Masset & M. Tagger: Propagation of warps in moderately thick disks 



This error is plotted on figure |[ We see that it is always 
very weak, so that the approximate dispersion relation is 
excellent to describe the "hidden" mode. This dispersion 
relation is a generalization of the one obtained by Nelson 
(1976) without shear and self-gravity, which was: 



v(r) 



cD^ 



4f}2 



+H/k 



For a ratio ag/a^ of about one fifth, we find a radial 
wavelength of about one kiloparsec, which might be varied 
by a factor 2 by varying the parameters. This is the order 
of magnitude of the wavelength of corrugations observed 
in our Galaxy and in some other edge on spirals (Quiroga 
et al. 1977, Florido et al. 1991). Since these corrugations 
are observed in the very young stellar population, i. e. are 
assumed to trace the motion of the gas disk, we consider 
that this "hidden" mode is a very good candidate to ex- 
plain them. 

One must note that the halo contribution to the ver- 
tical oscillation frequency should also be included in the 
approximate dispersion relation (see e.g. Toomre, 1983). 
This contribution depends only on the midplane density, 
so that for a given surface density of the halo it also de- 
pends on its scale height - of which nothing is known. Let 
us assume that the halo is a flattened disk in hydrostatic 
equilibrium. Let us also assume that its density dominates 
the local gravity, so that the stellar disk scale height is: 

Ph ii 
where the subscript H notes halo quantities, while 
an 

Then one easily finds that the Toomre parameters of the 
stars and the halo are in the ratio: 



c^^ 



m 



p* 

PH 



Thus, if both the halo and the stars have Toomre param- 
eters not too different from one, their midplane densities 
and their contributions to Vg must be comparable, so that 
our estimate of the corrugation wavelength remains valid - 
though this suffers the same uncertainties as any estimate 
concerning the halo. 



5. Amplification of m > 1 v^rarps 

In this section we briefly show that m > \ warps can be 
amplified, through a weak form of the Swing mechanism 
which amplifies spiral waves or modes. The m — 1 warp, 
which is of course the most interesting since it is the one 
observed in most galactic disks, is not concerned since 
it has been shown by Hunter and Toomre (1969) that it 
always has a positive energy, while spiral waves are ampli- 
fied by exchanging energy between negative energy waves 




Fig. 6. The propagation properties of bending waves are 
summarized in this diagram. As in the previous figures, v 
represents u> / k. It is a monotonically increasing function 
of r. See text for explanations. 



inside corotation and positive energy ones outside corota- 
tion. This is of course connected with the well-known fact 
that the shearing sheet model is not adapted to describe 
the m — \ mode, for which the effects of the cylindrical 
geometry cannot be neglected, (this is different from the 
mechanism studied by Bertin and Mark (1980), who in- 
vestigated the possibility of self-excited bending modes in 
galactic disks through the exchange of angular momentum 
between the disk and a slow bulge-halo component, with 
a "quantum condition" between the center of the disk and 
the corotation.) 

We have summarized on figure || the propagation prop- 
erties of bending waves. As for spiral waves, we have a 
forbidden band around corotation, imposed by the verti- 
cal resonances at -|-/x and — /i. In the case of spiral waves, 
the thickness of this forbidden band tends to zero as the 
Toomre parameter Q gets closer to unity. In the case of 
bending waves, this thickness tends to zero as the rota- 
tion curve tends to be flat (since p? = 2^^ — k^). The 
sign of the group velocity, which is the same as the sign of 
dv/dkx, allows one to give the direction of propagation of 
these bending waves. In our approximation (WKB and the 
perturbative derivation of the compressibility term), the 
latter term becomes infinite at the Lindblad resonances, 
so that we cannot conclude here on the exact properties 
of the warp close to these resonances. We reserve a more 
detailed study for future work, but we have checked that 
this does not affect the results presented here. 

Our motive for the investigation of non-WKB effects 
is that the dispersion relation of warps is similar in struc- 
ture to that of spiral waves, with a positive rather than 
negative self-gravity term. It is thus analogous to the dis- 
persion relation for spiral waves in a disk embedded in a 
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strong vertical magnetic field, as derived by Tagger et al. 
(1990). They showed that these waves were still subject 
to a weak form of the Swing amplification mechanism, as- 
sociated mathematically with the presence of the square 
root {q = (k'^ + kyY/^) in the self-gravity term. Physically, 
this corresponds to the long-range action of gravity or of 
magnetic stresses. Thus we can expect the same mecha- 
nism to apply to warps. 

However this amplification takes place when k^ is of 
the order of ky, so that the WKB approximation can- 
not be used straightforwardly. Although we might use the 
method of Pellat et al. (1990) for an analytical derivation, 
we present here a full numerical solution which allows us to 
derive the result directly. We stay in the shearing sheet and 
use the formalism developed by various authors (Goldreich 
and Lynden-Bell, 1965; Toomre, 1964; Drury, 1980; Lin 
and Thurstans, 1984), reviewed in Tagger et al. (1994), 
to compute the amplification. We refer to these works for 
a complete description of the formalism, where one uses 
the fact that because of the shearing motions kx changes 
with time as /c° — 2Akyt. Thus the uj term obtained in the 
WKB approximation, is replaced in the Eulcr equation by 
derivatives: 

This in turn can be, through a change of variables, trans- 
formed to a single derivative over k^. There results a set 
of differential equations which, in the tightly-wound limit 
{kx ^ ky), reduce to the WKB result while at low kx 
they give rise to transient phenomena and amplification. 
We write and solve numerically these equations. 

At large (positive or negative) kx one recovers the 



WKB results, while amplification can occur at kx 



fCi, 



The amplification can be described in the following man- 
ner: one can combine the equations and find at large kx 
WKB solutions, varying as 



M{kx) = P"^/*exp ±z / P^/'\k'Jdk'^ 



(13) 



where M is any perturbed quantity and P{kx) a quantity 
which appears in equations as the square of an "instanta- 
neous" frequency. These solutions correspond to leading 
and trailing waves respectively for kx negative or positive, 
propagating inside or beyond corotation respectively for 
the -|- and - signs in the exponent. The leading waves prop- 
agate toward corotation, and trailing waves away from it. 
Thus let us start at kx — > — oo (i.e. leading waves) 
with a "pure" solution (i.e. with the -I- sign in the expo- 
nent), corresponding to a wave traveling inside and toward 
corotation. The numerical integration over kx corresponds 
to following the wave as it approaches corotation, and then 
is reflected when kx becomes positive. But in the mean- 
time the WKB approximation has lost its validity while 
kx '^ ky, so that the "pure" solution has lost its identity. 



Thus in general, when the WKB approximation becomes 
valid again (at large and positive kx), we should obtain a 
mixture of the two solutions, with the -I- and — sign in the 
exponent: this means a mixture of trailing waves traveling 
inside and outside corotation, and away from it: as the 
initial leading wave is reflected it also "tunnels" through 
the corotation region, and generates an outgoing trailing 
wave. This is the usual process of the Swing amplifica- 
tion mechanism for spiral waves; in that case the complex 
physics of this mechanism implies that the reflected wave 
has a larger amplitude than the initial one, because waves 
have negative energy inside corotation and positive energy 
outside. Thus that by conservation of energy the emission 
of the positive energy wave outside corotation means an 
amplification of the negative energy one. 

In the case of warps the mechanism is exactly sim- 
ilar; however the repulsive, rather than attractive, force 
between density perturbations on cither side of corotation 
makes it much less efficient than in the case of spirals: one 
can obtain an amplification by a hundred in the most ex- 
treme case for spirals, while Tagger et al. (1990) find only 
a factor 1.05 (and a relative amplitude .3 for the transmit- 
ted wave) in the best case, when they study magnetized 
disks with a dispersion relation mathematically similar to 
that of warps. 

Here we obtain similar results, summarized in fig- 
ures {uh and (0). This ampliflcation is so weak that it 
is negligible for a single wave packet traveling in the disk. 
One might consider building a normal mode, formed as 
in the case of spiral when the reflected trailing wave is re- 
flected back at the galactic center to become a leading one 
traveling back toward corotation; then, with an ampliflca- 
tion by a factor 1.05, and assuming that the wave suffers 
no dissipation, the e-folding time would be about 20 times 
the duration of this cycle, which is at best a few rotation 
times. This means that such modes could not reach a siz- 
able amplitude in a galactic disk over a Hubble time. On 
the other hand, in accretion disks where the possible num- 
ber of rotation periods is much larger, and still neglecting 
all possible sources of dissipation, this might provide a 
way of maintaining long-lived warp modes. 

We wish to note that this result has also been obtained 
independently by J. Goodman (private communication). 

6. Discussion 

We have investigated the dispersion relation of warps in a 
moderately thick disk. We have found that the propaga- 
tion of warps is described by the dispersion relation (|lO|) , 
with the following hypothesis: 

— We work in the shearing sheet model; 

— The disk is consistent (with our deflnition, see the cor- 
responding section), which physically means that the 
disk must be moderately thick: it is geometrically thin 
but the ratio of the sound time through the disk to the 
rotation time is arbitrary; 
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Re[<Z>] 



F ig. Y . In this figure wc show the mean deviation of the warp from 
the galactic plane, i.e. I pQ(z)s{z)dz/ (12/2)^ solved from the full set 

of differential equations over k^. The parameters of the simulation are 
indicated in the diagram. We have set 47rG — 1 and H — 1. The units 
in Z are arbitrary, since the perturbation is linear. The vertical profile 
of p is imposed through the use of the consistency equation (9) by the 
choice of the parameters a, Pm and A/Q. Here, the high value of Pm (in 
our units) and the value of A/Q close to —1/2 give a disk dominated 
by the local gravity, i.e. with a vertical profile close to 1/ cosh (z/H). 
This can be easily checked with the values of the parameters /i/k and Q, 
which appear to be respectively, in this case, 7.1 X 10~ and 0.63 (just 
above the stability limit for this type of moderately thick disk). The solid 
line represents the real part, and the dotted one the imaginary part. For 
fca, < —ky, the solution remains WKB, with a decreasing frequency and 
an increasing amplitude (respectively due to the behavior of Cj with q and 
to the P~ ' term in the WKB solution). Since we have only one wave 
propagating in this region, real and imaginary parts are in quadrature. 
After crossing the region where fc^ ~ fcy, the real and imaginary parts 
are no longer in quadrature, because the solution now projects onto both 
WKB solutions, and we see that total amplitude has changed. 



— The warp wavelength is larger than the disk thickness; 

— The disk is isothermal; 

— the waves are WKB, i.e. tightly wound. 

The corrective term is due to the effect of compress- 
ibility, acting when the sound time through the disk is not 
much smaller than the rotation time. 

The associated horizontal motions are important be- 
cause they allow one to consider the possibility of a non- 
linear mechanism to continuously excite the warps: indeed 
the beat wave of two m ~ \ warps (or of a warp with itself) 
is an TO = 2 perturbation, whose potential and horizon- 
tal motion are even in z] thus they have the wavenumber 
and parity of a spiral. This means that a spiral can inter- 
act non-linearly with warps and excite them, in the same 
manner that it can do with other spirals. Non-linear cou- 
pling of bars and spirals has been studied by Tagger et al. , 
1987, and Sygnet et al. , 1988. This mechanism was found 
very efficient and explained the behavior of numerical sim- 



h ig. 8. In this figure we plot the imaginary part of the warp de- 
viation versus its real part, for the same simulation as figure (7). The 
"trajectory" first follows circles of slowly increasing radius, as the "pure" 
leading WKB solution travels toward corotation. When k^ has become 
positive, the trajectory becomes elliptic, because the solution is now a 
mixture of perturbations varying as exp ibi I P dk'^ . The measure of 
the ellipticity gives the ratio of the transmitted to the reflected wave. 



ulations of galactic spiral; in a forthcoming paper we will 
consider the possibility and the efficiency with which non- 
linear coupling of warps and spirals might feed warps from 
the energy and angular momentum carried by the spiral 
wave to the outer parts of the galactic disk. 

Another potentially important effect of the compres- 
sional term relates to the possible existence of warp modes, 
i.e. standing wave structures, in the disk. Classically, the 
Hunter- Toomre criterion shows that modes can exist only 
if the integral J q dr is finite, and that with the thin-disk 
dispersion relation this can be written as J dr/T, finite - 
a condition that is difficult to fulfill in realistic disk mod- 
els. With the new compressional term, this connection be- 
tween the two integrals disappears, so that one might have 
J q dr finite {i.e. a discrete mode spectrum) even when 
J dr/H is not bounded^. 

We have also investigated the self-amplification of 
m > 1 warps by a weak form of the SWING mechanism. 
We have found an amplification which is too weak to sig- 
nificantly amplify warps in galactic disks, but might be 
considered in accretion disks. 



^ As this paper was being rewritten following referee's com- 
ments, we learned that Sellwood (in preparation), using a mod- 
ification of the dispersion relation from Toomre (1966) has 
reached similar conclusions and indeed found standing warp 
modes in numerical simulations. 
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8. Appendix: derivation of the dispersion relation 

We derive hereafter the relation dispersion of warps in 
the monofluid case under the hypothesis mentioned in the 
main text. 
Let N be: 



after straightforward computations using the continuity 
equation (||): 



+ OC 



oz 



AttG 



+ 00 



PQadz- 



AnGq 



+ CXD 



Po{z)dz 



4-00 



adz + OgraviiqHy 



N = I a{z)dz 
'o 



+ OC 



po ( 2 ds 



uj \ dz 



dz 



(14) 



where use of equation dq) has been made. Using the 
unperturbed hydrostatic equilibrium and the expression 
of s given by equation (0) we rewrite TV as: 



u^ 



^^Sf^-d^da^^lJaYS 







dz dz 



-t-oo 



where we have defined 



ijj — K^ — a'^q 



d(j) 
oz 



(15) 



The a^q^ term in the denominator is the source of the 
new effect we introduce. It clearly represents the effect 
on s of the compressibility associated with the horizontal 
motions. 

Integrating by parts the first term, we use the Poisson 
equation and the consistency hypothesis (i.e. that A^^o 
does not depend on z), and get: 



where Ograv[{qH)^] means that this gravitational con- 
tribution also contains a term of second order in qH which 
we have not evaluated, truncating the expansion of the ex- 
ponentials to first order. An analytical computation of this 
can be performed if one assumes that the motion in the 
warp is vertical with a vanishing divergence, and if one 
has an analytical expression (e.g. a gaussian) for the ver- 
tical equilibrium density profile. In these conditions one 
finds that the dispersion relation is modified by replac- 
ing the surface density S by an apparent surface density 
S' = S(l — XqH), where A is a constant equal to {2n)~^^'^ 
in the gaussian case. Thus this term affects the dispersion 
'^f elation only in the limit qH ^ 1, whereas the compres- 
sional term we have derived can play a role for an arbi- 
trarily small value of qH. The computation of this grav- 
ity term could be well approximated by the use of "soft- 
ened gravity" models (Erikson 1974, Athanassoula 1984, 
Romeo 1994) where one artificially truncates the {r — r')~-^ 
divergence of the potential to take into account the geo- 
metric effect of the finite disk thickness. Thus hereafter 
we neglect this second-order gravitational term. 

Substituting our result and the value of S in equation 
([l6[), and dividing by N, we find the dispersion relation: 



LO^ 
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Where we have made use of the parity properties of 
the warp (s = and = at z = 0). 

We see here that the consistency hypothesis is the key 
to our derivation, since it has allowed us to extract Aj.0o 
from the integral, in the second term of the right-hand 
side. 

We write the solution of the Poisson equation (Equa- 
tion 7) as: 

- j^rGpoiz^Hz^l , 



-qz 



2q 

ATTGpa{z')s{z') 
2^ 



dz' 



3«^ dz' 



(one easily checks that this is the solution which veri- 
fies the boundary conditions ai z = and at z — > cxd). 

Since pq vanishes beyond a vertical scale H we can 
expand the exponential to first order in qH; this gives, 
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